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1 INTRODUCTION 



ABSTRACT 

A new method is presented to obtain a non-parametric maximum likelihood estimate 
of the luminosity function and the selection function of a flux limited redshift sur- 
vey. The method parameterizes the selection function as a series of step-wise power 
laws and allows possible evolution of the luminosity function. We also propose a new 
technique to estimate the rate of evolution of the luminosity function. This is based 
on a minimization of the observed large-scale power with respect to the evolutionary 
model. We use an ensemble of mock surveys extracted from a N-body simulation to 
verify the power of this method. 

We apply our estimators to the 1.2-Jy survey of IRAS galaxies. We find a far- 
infrared luminosity function in good agreement with previously published results and 
evidence for rather strong evolution. If the comoving number density of IRAS galaxies 
is assumed to scale cx (1 + z) p , we estimate P — 4.3 ± 1.4 ± 1. 

Key words: methods: statistical - galaxies: evolution - infrared: galaxies. 



likelihood estimates of the LF . In the so-called parametric 



maximum likelihood es timate (Sandage, Tammann & Yahil 



the larg^-sCalC structure of the Universe. However, sitch cat 
alogues show a strong decline of the mean number density of 
galaxies as a function of distance, simply because at larger 
distances more and more galaxies fall below the apparent 
flux limit. Most of the statistics that are used to analyse 
these surveys require an accurate knowledge of this depen- 
dence, which is described in terms of the selection function 
(SF). Hence the accuracy of the adopted SF can limit the 
reliability of such large-scale structure studies. 

Closely related to the SF is the luminosity function 
(LF), which describes the distribution in luminosity of the 
galaxy population sampled by the particular redshift survey. 
This quantity is of more fundamental importance, since it 



Flux limned redsmtt surveys are a major tool tor studying 197g |. [y ahi i et a l. 1 99 1|) an analytic form for the LF (or SF) 



is assumed that depends on a few (typically two to four) 
parameters. These are then determined by maximizing the 
likelihood of the observed data set. 

However, because the maximum likelihood technique of- 
fers no built-in measure of goodness-of-fit, almost any func- 
tional form can be made to 'fit', although the function may 
provide only a poor description of the data. The paramet- 
ric technique therefore requires an a priori knowledge of a 
suitable fitting form. 

If this information is not available one can allow a 
very flexible shape of the LF by describing it by many 
parameters in a reasonable way. This so-call ed step-wise 



should De reproduced by any viable theory of galaxy forma- 
tion and evolution. 



non-parametric maximum likelihood metho d (Nicoll & Se- 
Efstathiou, Ellis & Peterson 1988) has been used 



1983 



He 'e we revisit the problem of determining the LF and 



both to find luminosity functions (Efstathiou, Ellis fc Pe- 



terson 1988 



SF, given only the data of a flux limited redshift survey. Cur 



Loveday et al. 1992, Lin et al. 19961) and to 



rent standard methods for this task include Schmidt's (1968 ) 
VVKnax estimator, and the maxim um likelihood techniques 
first intr oduce d by Turner (L979) and Sandage, Tamman 
& Yahil (197E). The maximum likelihood methods are gen- 
erally superior to the older techniques because they allow 
the construction of estimators which are not systematically 
biased by density inhomogeneities. For this reason we will 
focus on them in the following. 

Two basic procedures have been used to find maximum 



estimate the run of radial densi ty with distance (paunder^ 



et al. 1990; Loveday et al. 1992). So far these applications 



only employed simple step functions to model the desired 
function. 

In this paper we propose a non-parametric maximum 
likelihood estimator that uses a new parameterization in 
terms of piece-wise power laws. This method provides ac- 
curate information on the shapes of the SF and the LF and 
has a number of computational advantages. For example, 
it does not require iterative solutions and it provides error- 
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estimates easily. The method is particularly useful for justi- 
fying specific analytic fitting forms for the SF and LF. 

We also discuss ways to estimate evolution of the LF. 
Such evolution appears quite strong in the far-infrared. A 
V/Vmax test may be us ed to estimate the evolutionary rate 
(Avni & Bahcall 198C), but the result can be stron gly af - 
fected by density inhomogeneities. Saunders et al. (199C) 



proposed an alternative maximum likelihood estimator for 
density evolution. Unfortunately, their technique is troubled 
by the same problem; it explicitly assumes a locally uniform 
distribution of galaxies. In order to reduce the influence of 
density inhomogeneities (which are clearly present) we pro- 
pose to combine the SF determination with a counts-in-cells 
analysis on large scales. The evolutionary rate is then esti- 
mated by minimizing the variance of the large scale density 
field. We demonstrate the robustness of this method using 
an ensemble of flux limited mock catalogues constructed by 
observing an N-body model universe. 

Our paper is organized as follows. Sections 2 and 3 
present in detail our non-parametric estimators for the SF 
and the LF. The estimation of evolution is discussed in sec- 
tion 4. Finally, we apply these methods in section 5 to the 
1.2-Jy survey of IRA S galaxies. In particular, we give results 
for the far-infrared LF, its evolution and the SF of the 1.2-Jy 
survey. 



2.2 The likelihood expression 

We now imagine that the catalogue is drawn from an under- 
lying parent distribution given by n z (r,L) for L > L m i n («) 
and by for L < L m j n (z), where z = z(r). Then the con- 
ditional probability p(Li\zi)dL that a source observed at 
redshift Zi falls into the luminosity range [Li, Li +dL] takes 
the form 



p{Li\zi) AL 



n Zi (ri,Li)AL 



n Zz (n,L')AL' 



(3) 



The denominator simply counts the available number of 
galaxies at that distance and the numerator gives the num- 
ber of galaxies in the particular luminosity range. Upon in- 
sertion of equations (|l]) and (|^) this becomes 



p(Li\zt 



$ H (Lj) 
S(z 4 ) ' 



(4) 



Note that the density fluctuations have dropped out of this 
expression and so the dependence on r can be dropped in 
equation Q. This insensitivity to density inhomogeneities 
makes the maximum likelihood technique used here supe- 
rior compared to older methods like the ordinary V/V max 
estimator. If one now maximizes the likelihood 



£ = TTp(Li[;Zi) 



(•») 



2 A NON-PARAMETRIC ESTIMATOR FOR 
THE SELECTION FUNCTION 

2.1 Definitions 

Let the field n z (r,L)AL describe the comoving number 
density of galaxies at epoch z, in the luminosity interval 
[L,L + AL] and at comoving spatial position r. Assuming 
that the luminosity distribution is independent of clustering 
the number density field may be written as 



« z (r,L) = ^ 2 (L), 



(1) 



where $z(L) describes the LF. Here n z (r) = n z (r, L) AL 
is the local number density and n z — (n z (r)) signifies the 
mean number density, averaged over many realizations of the 
Universe. The LF is normalized as f°° & Z (L) AL = n z and 
the dependence on z takes care of a possible time evolution, 
if present. The luminosity cut Lo may be used to handle a 
possible formal divergence of the integral over Q Z (L) at the 
lower end. 

We define the selection function S(z) as the mean co- 
moving number density of galaxies that one expects to see 
in a flux limited survey at redshift z. Then S(z) is given by 



S(z) 



$ Z (L)AL, 



(2) 



>0) 



where L m in(z) denotes the minimum luminosity a source at 
redshift z may have without falling below the flux limit of 
the catalogue. In this work we neglect the peculiar velocities 
of galaxies and take all redshifts to be cosmological, i.e. we 
adopt a simple redshift-distance relation z = z(\r\). 



of the whole data set with respect to S(z) one obtains an 
estimate of the SF that is not systematically biased by local 
density fluctuations. 

In order to find this maximum in practice we first ex- 
press $ Zi (Li) in terms of the SF with the help of equation 
(^) . Here a model for the evolution of the LF has to be spec- 
ified. For brevity we will only treat a case with pure density 
evolution according to 



*,(L) = g{z)$ {L). 



(6) 



However, our method can be easily generalized to include 
luminosity evolution as well. 

We further define a maximal redshift zf 1 for each source 
such that Lmin^f 1 ) = Li. If then the derivative S'(z) of 
equation (H) is evaluated at zf 1 one obtains 



s'(z?) = -iifl^L^LLnizT) + 4€r stf 



9(*r) 



(7) 



so that the probability of equation (|4j) can be expressed en- 
tirely in terms of S(z) and g(z). Hence one finally has to 
maximize 



A = ln£ = J^ln 



-S'{z?) + g'{zT) S{z?) 



S( Zi ) g{z?) S( Zi ) 



\g( Z f 



g(zi) 



InL^^D, (8) 



where the constant sum of the last term may be dropped. 

The above form suggests that one might be able to max- 
imize A simultaneously for g(z) and S(z). However, if equa- 
tion (Q) is rewritten for the density evolution model it be- 
comes 



p(Li\zi) = 



g(zi)$> {Li) 



$o(Li) 



AL 
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So the function g(z) drops out completely and there is no 
sensitivity to density evolution with this estimator. In fact, 
this was to be expected since estimates based on the like- 
lihood (Q) are independent of the density distribution by 
construction. 



2.3 A non-parametric maximum likelihood 
estimator 

Here we propose a new variant of the non-parametric 
method that models S(z) as a series of continuously linked 
power laws. This description seems appropriate since the SF 
is a smooth curve that covers a wide range of values and its 
local behaviour can be very well approximated by a power 
law. 

We describe S(z) by n pieces. Let S k = S(xk) be the 
values of S(z) at a series of ascending redshifts x k where 
k 6 {1, 2, . . . , n}. Then bin 1 covers < z < X\, bin 2 covers 
x\ < z < X2, and so forth. In each piece k, S(z) is taken to 
be a power law of the form 



S(z) = S k ( — 

s, Xfo 



for x fc _i < z < x k 



(9) 



where m k is the logarithmic slope of the particular piece. 
These slopes are precisely the quantities needed to charac- 
terize the shape of the SF. 

The different pieces have to join continuously, because 
S(z) is an integral over the LF. Continuity requires the Sk 
to be related by 



S k = Si J] 



J=2 



for 1 < k < n . 



(10) 



Let us further define at as the number of the bin to 
which the maximal redshift of galaxy i belongs. Similarly, 
let bi be the number of the interval that encloses the redshift 
Zi of galaxy i. Then the likelihood (H) takes the form 



A 



E 



In 



+ 



ff'K m ) 

<?or) 



+ 



In- 



E 



rribi In 



E 

j=6; + l 
Zj 
Xbi 



rrij In ■ 



Xj-l 



InLmin^i 11 )- 



(11) 



rrik, which might then quickly be improved by an iteration 
technique. This starting value can be calculated as 



n k g'jxk) 

Tk +Xk g{x k )' 



(14) 



where n k is the number of galaxies in bin k. Note that in the 
case of no evolution the solution to equation (|l^) is simply 
given by m k = -rikTu' 1 . 

The maximum likelihood method also allows an esti- 
mate oftlT£_£tatistiral_jrncertainties of the derived parame- 
ters (Kendall & Stuart 197E). Asymptotically the distribu- 
tion of C is a multivariate Gaussian around the true values 
of the parameters. If the information matrix / is defined as 



gHA 
lJ dmidrrij 



(15) 



then the covariance matrix of the parameter estimates is 
given by cov(mi,mj) = evaluated at the maximum 

of A. This result holds asymptotically for large sample sizes. 

Because the T k don't explicitly depend on any of the 
m k the information matrix is simply found to be 



0k, ai 



(rrik 



, g'(*n 



(16) 



One can therefore trivially invert this diagonal matrix and 
find error estimates for the m k as 



var(mfc) 



^. (17) 

n k 



This is a surprisingly simple result. In particular, one 
can solve for each of the m k independent of all the others. 
The m k are mutually uncorrelated which shows that it is 
essentially the local logarithmic slope of the SF that is de- 
termined by the maximum likelihood estimator ^j. 

Once the non-parametric estimate of the shape of the 
SF is found, it can be used to find and justify an appropri- 
ate analytic fitting function for S(z). For this purpose one 
can directly employ a minimum x 2 fit of the m k to the log- 
arithmic slope of some fitting form for S(z). Of course, once 
an analytic form has been selected in this way, the values 
of its parameters can also be determined with the paramet- 
ric maximum likelihood technique by using the fitting form 
directly in equation (|). 



The best estimates for the m k can be found by solving the 
likelihood equations 



OA 

dm k 



E 



m k 



+ T fc = 0, 



where T k is defined as 



(12) 



2.4 Normalization 

Because the normalization of the SF is lost with the above 
estimator it has to be found in a second step. For this pur- 
pose we write the SF as 



S(z)=i>s(z), 



(18) 



E E - 



Sj,k In ■ 



+ 



Cj-l 



ES k ,b, In — • 
x bi 



E Sk ' ai ln 



(13) 



The equations (|l^) are not fully linear in m k , but al- 
most. Since the redshift bins are narrow, replacing zf 1 by x ai 
will give a useful first approximation mj, to the true solution 



where s(z) is the shape of the SF as determined above. Then 
an unbiased estimate ip of the factor ip is given by 



f v m(r)w(r) dr 
J s(r')w(r') dr' 



(19) 



for an arbitrary weight function w(r). Here m(r) = 
y], 5(r — r{) represents the observed galaxy field and we 
employ the shorthand notation s(r) = s(z(|r|)). 
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Following Davis & Huchra (1982) we choose the weight 



function w(r) such that the expected variance of ip is mini- 
mized. This is to a good approximation the case for 



w(r) 



1 + J 3 ips(r) ' 



(20) 



where J3 = J 47rr 2 £(r) dr is the second moment of the two- 
point correlation function and V denotes the volume used 
in the normalization. 

Note that the normalization estimate has to be found 
iteratively since the estimator depends implicitly on ip. One 
also needs a model for the two-point correlation function in 
order to estimate J3. However, this is uncritical because the 
dependence on J3 is usually quite weak. Hence an estimate 
of Js based on a linear theory CDM power spectrum should 
be entirely sufficient. 

The numerical value of ip itself does not provide a mean- 
ingful measure for the local number density of galaxies. In 
fact, such a measure is not well defined for a purely flux lim- 
ited catalogue. Instead, the normalization is perhaps best 
expressed in terms of the expected number density N of 
galaxies per unit solid angle on the sky, i.e. 



N 



-L / s(A, 

47T / dz 



(21) 



3 A NON-PARAMETRIC ESTIMATOR FOR 
THE LUMINOSITY FUNCTION 

Using equation (^) it is possible to recover the underlying 
luminosity distribution if the SF is known, or vice versa. 
In particular one can readily derive a non-parametric LF 
estimate from the SF estimator described in the previous 
section. 

We suppose that the SF has been determined with the 
estimator described above and that the function L m i n (z) and 
its inverse z max (L) are known. The piecewise description of 
the SF directly translates into a piecewise description of the 
LF if boundaries Lk of luminosity intervals are defined by 
Lk = Lmi n (xk)- Upon evaluation of the derivative of equa- 
tion (0) at z — z max (L) the present day LF results as 



*o(£) = 



g'(z n 



- S (2 ma x ) S (2 ma x) 



g(z n 

X ^<7(£ m ax)L m i n (z m ax)) 

This translates into the piecewise description 



(22) 



$o(L) = 



g'{z n 



g(z n 



•2-rnax J V 3?fc 



(23) 



S , ( Z max)LJ nin (Zmax) 

for Lk-i < L < Lk, which is the desired non-parametric 
estimate. 

The usual way to find a non-parametric maximum like- 
lihood estimate of the L F utilizes a parameterization by a 



series of step functio ns ( Efstathiou, Ellis fc Peterson 1988 
1992| ) 



Loveday et al. 1992Q . We see two main advantages of our 
method as compared to this approach. 

First, the shape of the LF over each bin is approx- 
imately a power law, so it is able to adapt to the true 



shape of the LF in a flexible way. Only very small dis- 
continuities at the boundaries of bins remain and the es- 
timate &o(L) of equation ( |23[ ) traces the smooth LF quite 
well. The conventional estimator on the other hand assumes 
constant $0 over each bin which leads to large discontin- 
uous jumps in the LF at bin boundaries. Efstathiou et al. 
( 1988| ) show that the heights <3>fc of these bins are related 
^ ^ « K dN(L) to the un- 

derlying luminosity distribution. It is therefore more com- 
plicated to infer the smooth underlying LF from this non- 
parametric estimate. 

A second nice feature is that the absence of correlations 
among the rrik simplifies the estimation of errors in the LF. 
To demonstrate this let us define the values = &o(Lk) 
of the LF at the points Lk- In appendix A we compute the 
covariance matrix 



cov(]n$ fc ,ln$i) 



(24) 



of these quantities. The matrix Vki may then be used to fit 
an analytic model to the non-parametric LF estimate by 
minimizing the covariance form 



X 2 = ^(ln$ ; 



•lu#,)^» (In** -In**). 



(25) 



Although not of major importance a further practical 
advantage of the estimator presented here is its relative com- 
putational ease. In particular one does not have to resort to 
lengthy iteration techniques as required by the conventional 
parameterization with a series of step functions. 



4 ESTIMATING EVOLUTION 

4.1 The radial density distribution 

As was demonstrated above, the estimators based on equa- 
tion (Q) are insensitive to density evolution. If the latter is 
important (as appears to be the case in the IRAS surveys) 
a different approach is needed to determine its rate. 

Instead of p(Li\ri) as before we can equally well write 
down the conditional probability 

dr = ^fr.KLQdr = 47rr?n, i (rQdr 
ft* nZi (r,Li)dV /;• n z (r)dF 

to find a galaxy with luminosity Li in a distance interval 
of width dr at coordinate r< < fl max . Here the density 
field has been averaged over direction and the definition 
r* = min(r™, 7i max ) allows an upper distance cut-off 7i ma x 
to be included. We consider the inclusion of the latter to 
be important because the density field is only known well 
to a finite depth. At large distances it is dominated by shot 
noise due to increasingly sparse sampling. In addition the 
completeness of most surveys is worst in this regime. 

Equation ( pj| ) shows that the density distribution n z (r) 
can be estimated independent of the LF if 



A = \np(n\Li) = ^ In 



/</ n z(r) (r)dV 



(27) 



is maximized with respect to n z (r). Saunders et al. (199C) 
have used a parameterization of n z (r) in terms of a series 
of step functions leading to a non-parametric radial density 
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estimator. This estimator is quite useful since it can indicate 
the presence of evolution and it can show the most promi- 
nent density features. 



4.2 The constant density method 

In order to ob tain a n estimate of the evolutionary rate 
Saunders et al. ( 1990 ) proposed replacing n z (r) in equation 
(^) by the expected mean density n z = nog (z) at epoch 
z = z(r). Adopting an evolutionary model, for example 



g(z) = (l + z) p , 

one can then find the rate by maximizing 
(l + z.f 



A 



n (1 



(28) 



(29) 



with respect to P. A number of subsequent studies (Fisher 



et al. 1992 



Oliver et al. 1994) also applied this method to 
IRAS galaxies. 

Of course, in reality the galaxies are not distributed uni- 
formly. This fact leads to a serious drawback of the above 
estimator; the resulting estimate of P can be strongly influ- 
enced by density inhomogeneities. For example, an overden- 
sity in the foreground mimics negative evolution, whereas an 
overdensity in the background can bias the evolution low. 

Since this effect should be strongest nearby, a lower red- 
shift cut-off z m in, as introduced above, may be used to partly 
reduce its influence. However, as will be seen in our appli- 
cation to the 1.2-Jy survey the evolutionary estimates can 
depend quite strongly on the particular redshift interval cho- 
sen and it is unclear which choice gives the most reliable 
answer. 



4.3 The minimal variance method 

We here want to propose an alternative evolutionary estima- 
tor that avoids the unphysical neglect of density fluctuations 
and that gets rid of the arbitrary choice of redshift interval 
used in equation (]29|). 

As shown above the data of the redshift survey alone do 
not allow us to disentangle density evolution and large-scale 
structure. This is only possible if the fundamental assump- 
tion of statistical homogeneity is added. In a pragmatic ap- 
proach we might then base the determination of evolution 
on the notion that the correct rate should lead to a universe 
that looks as close as possible to homogeneous on very large 
scales. As a quantitative measure for deviations from ho- 
mogeneity one might, for example, take the variance a 2 of 
the density field smoothed on some large scale A. The deter- 
mination of the evolution might therefore be formulated in 
terms of an extremal principle: The evolutionary rate can be 
estimated by minimizing the observed power on large scales. 

To estimate a 2 we apply a slight ly modified version of 



the method of Saunders et al. ( 1991 ). We first compute a 
smoothed galaxy density field 



di = > Wi 



i 



S, 



(30) 



smoothing kernel, which we here take to be a Gaussian of 
the form 



(31) 



Then unbiased estimators for the mean d and the vari- 
ance a 2 of the galaxy density field are given by 



d = 



and 



(2) 



where 



(n) 



(32) 



(33) 



(34) 



These estimators become minimum variance estimators in a 
good approximation if the weights are chosen as 

-i 



9i = 
and 

hi 



(2) 



(35) 



yw+2(y«)' 



+ 2a" 



3(K«V + 4Y/ 3 >+4Y/ 2 ^- 2 



(36) 



Since the estimate of the variance depends on the SF it 
becomes a function of the adopted evolutionary model used 
in computing the SF. It is therefore possible to estimate the 
evolutionary rate by minimizing a 2 . 

A difficulty with this technique is that it does not easily 
supply an error estimate for the evolutionary rate. Addition- 
ally it is unclear how to choose A in an optimal way. A large 
A allows the survey volume to be probed to greater depth, 
thereby increasing the sensitivity to evolution. However, the 
number of independent smoothing volumes declines simul- 
taneously, thereby degrading the accuracy of the variance 
estimate. What is then the optimal choice? 

We solve both of these problems by extracting an en- 
semble of mock catalogues from a large N-body simulation. 
The mock surveys mimic the statistical properties of the 
survey at hand in terms of the SF, the evolution, the sky 
coverage and the source density. This suite of catalogues can 
be used to find the optimal A and to estimate the precision 
of the final evolutionary estimate. 

Although computationally somewhat lengthy this 
scheme is well worth the effort because it offers an increased 
accuracy as we will demonstrate with the mock surveys in 
our application to the 1.2-Jy survey. A further advantage 
is that it is equally applicable to luminosity evolution, in 
contrast to the estimator discussed in section 1.2. 



APPLICATION TO THE 1.2-JY REDSHIFT 
SURVEY 



on a fine mesh. Here rrij denotes the number of galaxies in The data of the 1.2-Jy redshift survey (ptrauss et al. 199C) 



cell j, Sj is the value of the SF and Wi 



of IRAS galaxies has been published (|Fisher et al. 199 



I*) 
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and can be r etrieved electronically f rom the Astronomical 
Data Center (ftp://adc.gsfc.nasa.go\). The 5321 galaxies of 



the survey are selected from the PSC catalogue above a flux 
limit of 1.2 Jy in the 60 fim band. The sky coverage is 87.6 
per cent, excluding only the zone of avoidance for |6| < 5° 
and a few unobserved or contaminated patches at higher 
latitude. 

We convert the redshifts to the Local Group frame and 
use them subsequently to infer distances without further 
corrections for peculiar velocities. 



Table 1. Parameters of the selection function fit. 



0.741 



+0.128 



4.210 



-0.419 
-0.344 



1.582 



7 

+0.237 



0.0184 



z 

+0.00213 



if> [h 3 Mpc" 3 ] 
(486.5 ± 13.0) x 10" 6 



5.1 K-correction and maximal redshift 

The estimators presented above demand an accurate compu- 
tation of the maximal redshift zf 1 — z max (Li) that a given 
source could have without falling below the flux limit. A 
detailed calculation of this quantity involves K- and colour- 
corrections, which require a knowledge of the spectral en- 
ergy distribution (SED) of the source, a specification of the 
responsivity R(v) of the detector, and a choice for the cos- 
mological background model. 



For the case of the IRAS surveys Fisher et al. (1995) 
have phrased this problem in terms of the useful equation 



/mill 



r? 



V 



(37) 



which implicitly specifies zf 1 in a general form. Here /, is 
the quoted flux density of a source and rt is its comoving 
coordinate distance. / m i n gives the flux limit of the catalogue 
and r\ is defined as rj = (1 + zf)(l + The function 



J R(v)f v {vrj) Av 
fR(v>).U(v')dv> 



(38) 



encodes K- and colour-corrections by means of the observed 
flux density f v (v) of the source. The cosmological back- 
ground model enters via the function r = r(z). We will adopt 
an Einstein-de-Sitter universe, i.e. 



r{z) = 



ao Ho 



1 - 



1 



(39) 



The 60 /im band SED of IRAS galaxies might be ap- 
proximated as a straight power law L v (u) oc v a with 
a « —2. In this case *&(??) = rj a . We have typically used 
q = —2, but also tried a gray-body model (Saunders et al 



199C ) fitted to 60 ^m/100um flux ratios and the polynomial 



model of Fisher et al. ( 1992 ) which results in a considerably 
shallower SED. 



5.2 Selection function 

In our determination of the SF and LF we will assume that 
the LF exhibits density evolution of the form g(z) = (l + z) p 
with P = 4.3. The actual esti mat ion of the evolutionary rate 
is discussed below in section 5.4. 

Figure |l| shows the estimated local slopes of the SF 
and a fit based on the functional form 



S{z) 



(i + (^) 7 ) 



(40) 



We used 40 redshift intervals, logarithmically spaced be- 
tween xq = 0.003 and X39 = 0.15, for the non-parametric 




0.05 0.10 
redshift z 




0.05 0.10 
redshift z 



Figure 1. Non-parametric selection function estimate. The top 
panel shows the estimated slopes with la error bars. The fit 
is based on equation (fi"o|). The lower panel compares the actual 
non-parametric SF (solid) and the analytic fit (dotted) with the 
parameters of table 1 . 



estimate. A minimum \ 2 fit of the measured mk to the form 
of equation ( |40| ) resulted in a reduced xt °f 0-87 (for v = 36 
degrees of freedom), which indicates a good fit. This form of 
the S F is sl ightly more genera l than the one used by Yahil 
et al. ( |l99l| ) and Fisher et al. (|l99§) who essentially fixed 7 
at the value 2. This gives a marginally worse \v of 0.91. 

Our best estimates for the final SF parameters are listed 
in table They are obtained by maximizing the likelihood 
of the parametric form directly. Although a fit to the non- 
parametric estimate gives a very close result, we prefer these 
numbers as final estimates, because they are free of any bin- 
ning effects. For each of the parameters a, /3, 7, and z* the 
quoted errors give la confidence intervals obtained from the 
bounding box of the A max — 0.5 likelihood cont our. 

Figure shows the ratio of Fisher et al.'s (1995) SF to 



our estimate. The agreement is quite good out to redshift 
around 0.1. However, relative to our result the predicted 
mean density of galaxies at re dshift 0.15 is about 20 per 
cent higher in the Fisher et al. ( [t995| ) result. 

In order to normalize the SF we used the volume inside 
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0.08 
red shi ft z 



Figure 2. Ratio of the 1.2-Jy-SF of Fisher et al. ( JL995[ ) to the 
result of this work. The ratio is scaled such that both parame- 
tcrizations of the SF predict the same number of galaxies out to 
redshift 0.15. 



z = 0.15. This determines the parameter ijj and results in 
TV = (490 ± 13) sr _1 expected galaxies per unit solid angle 
on the sky. Our final LF is also normalized to this num- 
ber. The la error-estimate is based on a direct computation 
of the variance of ip by means of equation (|ig|). Here the 
uncertainty in the shape of s(z) is neglected. 



5.3 Luminosity function 

Figure |^ displays our non-parametric estimate of the far- 
infrared LF resulting from the method outlined in section 3. 
For comparison al so sh own is the non-parametric estimate 
of Saunders et al. (199C ) which takes the form of a histogram 
because it uses step functions to model the LF. In order to 
facilitate a comparison with the literature we present the LF 
in terms of luminosity per decimal decade of luminosity, viz. 



P 10- 4 - 



4>{L) = $ (L)L In 10, 



(41) 



and we define L for a source with SED L v (v) as L = v L„ (y) 
at 60 fim. 

We believe that our parameterization offers an improved 
description of the LF, without large unphysical jumps due 
to binning of sparse data. In particular, the quality of a fit 
with an analytic function can be judged quite easily. 

As a fitting form we use the function 



L 



(42) 



proposed by Lawrence et al. ( 198f ) and give the best-fit 
parameters in table |j| The fit is based on a minimization 
of the covariance form (|2^), imposing a fixed normalization 
as described in appendix A. The cited errors for a, j3 and 
L t give 68.3 per cent confidence intervals derived from an 
increase of A^ 2 = X 2- Xmin by 1, where A\ 2 is marginalized 
in the space of a, f3, L*, and C. Figure |i| shows contour levels 
of constant Ax 2 in the subspace of a and /3. The contour of 
Ax 2 = 2.3 encloses the 68.3 per cent joint confidence region 
of q and (3. The error given in table ti for the normalization 
C is computed for a fixed shape of the LF. 



L [L a h1 



Figure 3. Luminosity function estimates. The top panel com- 
pares our non-parametri c LF estimate (thick) with the esti- 
mate of Saunders et al. ( L99fj| ). The latter has been rescaled to 
Hq = 100/ikms -1 and shifted vertically for graphical clarity. The 
two upper curves (also shifted) compare the non-parametric esti- 
mate with our analytic fit. 

The lower panel shows various parameterizati on o f the far- 
infrared LF by di fferent authors: Saunders et al. jl99Cj ) (dotted), 
Yahil et al. (1991) (dot-dashed), Lawrence et al. (119861) (dashed), 



this work (thick). 




1.05 1.10 1.15 1.20 1.25 1.30 1.35 



Figure 4. Contour levels of constant A\ 2 for the LF fit in the 
subspace of a and ft. The thick contour defined by Ax 2 = 2.3 
encloses the 68.3 per cent joint confidence region of a and j3. 
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Table 2. Parameters of the luminosity function fit. 

a 



1.221 



^0.068 
-0.072 



-+0.640 



9 1 lfi+0.080 



C [/i 3 Mpc" 3 ] 

(1.670 ±0.045) x 10" 2 



The simple-two power law Q42|) does a remarkably good 
job in fitting the measurements as evidenced by the compar- 
ison in figure ^ and by a reduced xt °f 1-06. In particular we 
find no need to choose a functional form that shows more 



curv atu re over the full range of luminosities ( Saunders et al i 
199C |). f lowever, we must not forget that we have ignored 
velocity fields in t his pa per and that we cannot go as faint 
as Saunders et al. fll990p who used an additional set of very 
local galaxies. This might affect the faint end slo pe; w e find 
a value not quite as shallow as Sa under s et al. ( |1990 ), but 
somewhat flatter than Yahil et al. ( 1991 ). 



5.4 Evolution 

Evidence for evolution in IRAS galaxies has been reported 



(Hacking & Soifer 19931; 


Hacking, Condon & Houck 1987 


Lonsdale et al. 19901; Gergorich et al. 1995) and redshift 


surveys (Lonsdale & Hacking 1989]; 


Saunders et al. 1990 


Fisher et al. 1992; Oliver et al. 1994 


). However, there has 



been some controversy about the magnitude of the evolu- 
tionary rate seen in the far-infrared LF. Saunders et al. 
( [i990| ) foun d P = 7 ± 2 for the QCD survey, whereas Fisher 
ct al. (1992) gave P — 2 ± 3 for the 1.2-Jy survey, a resul t 



apparently consistent with no evolution. Oliver et al. (1994) 
found P — 6.25 ± 1.5 for a survey of faint IRAS galaxies and 
P = 4 2 ± 2 .3 for the QDOT survey. More recently Bertin 
et al. fll996j ) report P = 6.0 ± 1.2 for the FIR sample. 



The low resul t for the 1.2-Jy survey, obtained with the 
was based on an early version of the 



method of section 4.2 



catalogue. Since evolutionary estimates are quite sensitive to 
completeness we here analyse the final catalogue again. This 
is also necessary because the evolutionary rate is needed as 
input parameter for the determination of the SF and LF. 



5.4-1 The run of radial density 

A first indication that a correct treatment of evolution is 
important may be obtained with the radial density esti- 
mator based on equation (|2^). We adopt a parameteriza- 
tion of n z (r) in terms of 45 step-functions (with i? max = 
400 /i _1 Mpc) of widths chosen such that the signal-to-noise 
ratio in each bin stays roughly constant. 

In the resulting distribution, displayed in figure ^, the 
average density seems to rise with distance. A similar be- 
haviour is also found if smaller patches of the sky are exam- 
ined separately. As an illustration we have split the sample 
in 4 regions of approximately equal size selected in terms of 
the modulus |6| of the galactic latitude. In this way we also 
obtain a low latitude sample which might be less complete 
than the other parts of the survey. In detail these regions, 
labeled B1-B4, cover |6| < 18° (Bl), 18° < |6| < 33° (B2), 
33° < |6| < 52° (B3), and 52° < |6| (B4). As evidenced 




200 

distance [Mpc/h] 






B2 : 




18°<|b|<33 = 



100 1 
distance [Mpc/h] 



distance [Mpc/h] 





100 1 
distance [Mpc/h] 



52 a <|b|<90° 

150 200 



distance [Mpc/h] 



Figure 5. Radial density distribution obtained with a non- 
parametric estimator based on equation (j2?j) and the assumption 
of no evolution. The result for the full sample is shown in the top 
panel, whereas the plots labeled B1-B4 are computed for subsam- 
ples selected by galactic latitude: |6| < 18° (Bl), 18° < |6| < 33° 
(B2), 33° < |6| < 52° (B3), and 52° < |6| (B4). 



Table 3. Evolutionary estimates P obtained with the constant 
density method for various redshift intervals. The subsamples Bl- 
B4 are selected in terms of galactic latitude: \b\ < 18° (Bl), 18° < 

161 



< 33° 


(B2), 33° < |6| < 52° (B3), 


and 52° < |6| 


(B4). 




[0,oo] 


[0.01, oo] 


[0,0.1] 


[0.01,0.1] 


Full 


3.7 ±1.5 


5.0 ± 1.6 


3.9 ± 1.9 


5.9 ±2.0 


Bl 


-1.0 ±3.3 


-1.5 ±3.4 


2.6 ±3.9 


2.5 ±4.2 


B2 


9.2 ± 3.2 


6.1 ±3.2 


11.1 ± 3.8 


6.7± 3.9 


B3 


6.1 ± 3.1 


7.8 ±3.3 


7.1 ± 3.7 


9.7 ±4.0 


B4 


0.6 ± 2.8 


6.5 ±3.1 


-6.2 ±3.8 


4.1 ±4.3 



by figure |H| the rise of density with distance can be seen in 
all four subsamples, albeit with variable strength. Promi- 
nent structures are also visible, in particular Perseus-Pisces 
(around 50/i _1 Mpc) and the Local Supercluster, which is 
responsible for the nearby overdensity seen in B4. 



5.4-2 Evolutionary estimates with the constant density 
method 

We now turn to estimates of the density e volu tion parameter 
P obtained with the estimator of section iJ2 . Table ^ shows 
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results for various redshift intervals, in each case computed 
for the full catalogue and for the subsamples B1-B4. 

Although the statistical errors are large, they are in- 
sufficient to explain the large differences in the estimates 
that result when the redshift boundaries are changed. For 
example, when the upper redshift cut-off is lowered from oo 
to O.f, the low latitude sample Bl shows a marked increase 
in P which is likely to be the result of incompleteness at 
large z. The sample B4 on the other hand exhibits a very 
strong overdensity nearby, that biases evolution low if it is 
not excluded by imposing a lower redshift cut-off. 

We have tried many different redshift intervals and find 
it hard to come to definite conclusions with this method, 
simply because the result is strongly affected by density in- 
homogeneities and the choice of a redshift interval exhibits 
a certain degree of arbitrariness. We therefore prefer to esti- 
mate evolution with the minimum variance method outlined 



in section 4.5 



5-4-3 An ensemble of mock surveys 

In order to demonstrate that this method is indeed more re- 
liable we have constructed an ensemble of f .2-Jy-like mock 
catalogues by observing a large N-body simulation of a stan- 
dard CDM universe. The simulati on has been kindly pro- 
vided by the Virgo collaboration ( Jenkins et al. f996 ). ft 
contains 256 3 particles in a periodic box of size V = L 3 
with L — 240/i _1 Mpc and is normalized to fluctuations of 
as — 0.6 in spheres of 8/i _1 Mpc. The details of this model 
are not of much relevance for the following, the salient point 
is that the model exhibits clustering in realistic strength. 

By randomly placing an observer in the periodic box we 
constructed 50 mock catalogues of this simulation. We used 
a SF very similar to the 1.2-Jy survey (with parameters a — 
0.84, (3 = 3.96, 7 = f .74, z* = 0.018) and imposed density 
evolution with P = 5. The mock catalogues featured the 
same angular mask as the 1.2-Jy survey and were sampled 
to a density of 532f galaxies inside 460 h~ Mpc. The suite of 
these catalogues simulates Poisson sampling noise and - to a 
reasonable degree - cosmic variance as well. The catalogues 
can be ideally used to test the performance of the various 
estimators discussed in this paper. 

Firs t we applied the constant density estimator of sec- 
tion 4.2 to these catalogues. Using the redshift interval 
[0,0.1] we obtained (P) — 5.7 with an rms scatter of 2.4 
among the individual measurements. This scatter is some- 
what larger than the statistical error of ~ 1.9 provided by 
the likelihood method itself. Smaller redshift intervals also 
led to estimates (P) close to 5, however with larger scatter. 

Next we applied our new minimum variance method to 
the mock catalogues. For this purpose we computed max- 
imum likelihood SF fits for 11 different values of P in the 
range to 10 individually for each catalogue. This allows 
the computation of a curve (t\(P) for each catalogue, and - 
by fitting the minimum with a parabola - an estimate of P. 

We treated the technical difficulty of the presence of an 
angular mask in the surv eys by using the ratio method of 
Melott & Dominik (1993) in the smoothing process. This 
means that we actually chose 



as kernel in equation (p(]|), where M(r) is a field equal to 1 
in the actual survey volume and equal to in the volume 
behind the angular mask. The estimation of the variance is 
only done with cells inside the survey volume. 

The method finally requires the choice of a smoothing 
length A. Since the minimum variance estimator automat- 
ically turns off when the signal-to-noise ratio becomes too 
low a larger A usually means that the the survey is probed to 
higher depth, thereby increasing the sensitivity to evolution. 
On the other hand, the accuracy of the measurement of the 
variance degrades with increasing A because the number of 
independent smoothing volumes in the usable survey volume 
declines. We find a choice for A that represents a compromise 
between these two effects by minimizing the rms scatter of 
the P measurements for the mock catalogues. 

For this purpose we have computed estimates of (P) and 
its scatter for a range of smoothing lengths and find that the 
uncertainty is smallest for A ~ 60 h~ Mpc, however, any A 
in the range 50-80 h~ Mpc works almost equally well. For 
A = 60ft _1 Mpc we obtain (P) — 5.3 with rms fluctuations 
of 1.4. This is a considerable improvement compared to the 
constant density method, nearly cutting the random error 
in half. This is possible because the new method eliminates 
most of the systematic influence of density inhomogeneities 
on the estimate of P. The remaining uncertainty is mainly 
due to counting statistics. 



5-4-4 Evolutionary estimate with the minimum variance 
method 

The above results show that the new estimator for density 
evolution gives a more precise estimate than the constant 
density method. Hence we apply it now to obtain our best 
estimate for the density evolutionary rate of the 1.2-Jy sur- 
vey. Adopting A = 60 /i _1 Mpc and using the same procedure 
as applied to the mock surveys we find P — 4.31± 1.4, where 
the error estimate is based on the mock surveys. 

Clearly, this error estimate is somewhat too optimistic, 
because the 1.2-Jy survey is not as perfectly selected as the 
mock catalogues. Additionally there are a number of sys- 
tematic uncertainties that degrade the accuracy of the evo- 
lutionary measurement. 

For example, the SED might be shallower as we assumed 
here, leading to an overestimate of evolution. The gray-body 
model (we find P — 3.9) does only very little change com- 
pared to the straight a — —2 SED used here. On the other 



hand, Fisher et al.'s (1992) significantly flatter polynomial 
model lowers the evolutionary rate quite strongly to P = 3.0. 
However, the SED of this model is likely to be too flat. In 
addition to 60 /im and 100 /im it tries to simultaneously fit 
the 25 fim fluxes which are typically underpredicted by the 
gray-body model. However, only half of the galaxies (the 
bright ones at 25 /jm) have detections in this band at all, 
so the polynomial model is biased towards a shallow SED. 
The polynomial model will also underpredict the slope of 
the SED if the spectra are indeed comprised of a cool and 



a warm component, as seem s often to be the case (Rowan- 
Robinson fc Crawford 1989[). 



f W{n-r>)M(r')Ar' 



(43) 



Because the 1.2-Jy survey is quite local, a change of 
background cosmology to a low density universe has very 
little effect, increasing the estimated P very slightly. 

More serious are potential problems with the IRAS 
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flux scale. Either a baseline flux error with accompanying 
Malmquist bias, or a nonlinear flux scale could lead to a sig- 
nificant overestimate of the evolutionary rate. However, so 
far there is no convincing evidence that the IRAS PSC is 
troubled with flux errors of the strength required to explain 
a significant fracti on of the evolutionary signal. For exam- 
ple, Oliver et al. (1994) give an upper limit of 36mJy for 



the baseline flux error of the PSC. At this level, Malmquist 
bias is not an issue, as we have checked with mock surveys 
that exhibit artificial flux errors. 

We believe that our new method for estimating the 
evolutionary rate has removed most of the uncertainty due 
to density inhomogeneities. However, sample incompleteness 
can bias the evolution low if a fraction of the faint high red- 
shift objects is missed. At low latitude, this is indeed quite 
likely the case for the 1.2-Jy survey. 

We think that an error of ±1 represents a reasonable 
estimate of these systematic uncertainties. Our final esti- 
mate for evolution is therefore P = 4.3 ± 1.4 ± 1. Given the 
uncertainties, this result is in good agreement with other de- 
terminations of P cit ed ab ove. However, is is notably higher 
than Fisher et al.'s (1992) result for an early version of the 
1.2-Jy survey. 



6 DISCUSSION 

An accurate determination of the SF of a redshift survey 
is a prerequisite for taking full advantage of its information 
about the large-scale structure of the Universe. For example, 
statistical methods that rely on estimates of the density field 
in large survey volumes, like power spectrum measurements 
or genus statistics, depend crucially on a precise knowledge 
of the mean density as a function of distance. 

In this work we have proposed a flexible non-parametric 
maximum likelihood estimator for the SF and LF. The 
method is independent of density inhomogeneities, gives ac- 
curate information on the shapes of SF and LF, and provides 
estimates of the statistical uncertainties of the derived quan- 
tities with relative computational ease. We think that the 
technique should be useful for upcoming redshift surveys. 

We have also proposed a new method to estimate evo- 
lution of the LF, based on the notion that the Universe 
should look homogeneous on large scales. With an ensemble 
of mock surveys we have demonstrated that the uncertain- 
ties due to density inhomogeneities, which troubled previous 
estimators, can be greatly reduced in this way. 

In our application of these estimators to the 1.2-Jy sur- 
vey of IRAS galaxies we have found evidence for strong evo- 
lution, confirming reports by several previous authors. Ex- 
pressed in terms of density evolution oc (1 + z) p , we find 
P — 4.3 ± 1.4 ± 1. This high evolutionary rate for the far- 
infrared LF is hard to explain as a statistical fluke. If con- 
firmed the strong evolution of the far-infrared LF will repre- 
sent an interesting challenge for theories of galaxy formation, 
interaction and evolution. 
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APPENDIX A: COVARIANCE MATRIX OF 
THE NON-PARAMETRIC LUMINOSITY 
FUNCTION ESTIMATE 

In this section we compute the covariance matrix 

V k i = cov(ln$ fe ,ln<I> i ) (Al) 

of the non-parametric LF estimate. Because we are primar- 
ily interested in the uncertainty of the shape of the LF we 
impose a fixed normalization, i.e. the quantity 



N' 



S(z)z 2 dz = SiQ 



(A2) 



is held constant. Here Q = Q(mi, . . . , m n ) is given by 

, n k 



E 

k=l 



rn k + 3 



-fir) n 

j=2 



Xj-1 



(A3) 
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Then ln$fc may be expressed as 



ln$t 



hi 



g'(x k ) m k 



In Q + In 



+ > m i m — — 



3=2 



N' 



g(xk)L' m ^(x k y 



(A4) 



An expansion of In to linear order in the slopes rrii allows 
an estimate of the covariance matrix in the form 



Vm = ^ A ki A u v&r(mi), 
i=i 

where we have defined 
91n$ fc 



drrii 



(A5) 



(A6) 



The expressions for A ki are somewhat lengthy, yet straight- 
forward to calculate: 



A k i 



S k i 



m k - Xk 



77 — r + ^ Sij In 



(A7) 



_L_ 

Q 



Ui- 



Si 



Si (m, + 3) 2 



Here f/j is given by 



E, i = In _5i_\-|L_5L_ 

Xi-i Si mi + 3 



for i > 1 and by Ui = for i = 1. 



m; +3 



(A8) 
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